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Abstract 

Some epidemics have been empirically observed to exhibit outbreaks of all possible sizes, i.e., to be scale- 
free or scale-invariant. Different explanations for this finding have been put forward; among them there 
is a model for "accidental pathogens" which leads to power-law distributed outbreaks without apparent 
need of parameter fine tuning. This model has been claimed to be related to self-organized criticality, and 
its critical properties have been conjectured to be related to directed percolation. Instead, we show that 
this is a (quasi) neutral model, analogous to those used in Population Genetics and Ecology, with the same 
critical behavior as the voter-model, i.e. the theory of accidental pathogens is a (quasi)-neutral theory. 
This analogy allows us to explain all the system phenomenology, including generic scale invariance and 
the associated scaling exponents, in a parsimonious and simple way. 

Introduction 

Many natural phenomena such as earthquakes, solar flares, avalanches of vortices in type II supercon- 
ductors, or rainfall, to name but a few, are characterized by outbursts of activity. These are typically 
distributed as power-laws of their size, without any apparent need for fine tuning ~ i.e. they are gener- 
ically scale invariant- [THB]- This is in contrast to what occurs in standard criticality, where a control 
parameter needs to be carefully tuned to observe scale invariance. The concept of self- organized critical- 
ity, which generated a lot of excitement and many applications in different fields, was proposed to account 
for generic scale invariance, i.e. to explain the "ability" of some systems to self-tune to the neighborhood 
of a critical point [lH3] . 

The spreading of some epidemics, such as meningitis in human populations, has been repeatedly 
reported to exhibit scale-invariant traits, including a wide variability of both durations and sizes of 
outbreaks. Moreover, the ratio of the variance to the mean of the distribution of meningitis and measles 
outbreak sizes have been empirically found to be very large and to grow rapidly with population-size [7] . 
This is the hallmark of anomalously large fluctuations such as those characteristic of the heavy tails 
of power-law distributions [5]. Actually, power- laws have been proposed to fit the statistics of some 
epidemics such as measles, pertussis or mumps in some specific locations as the Faroe islands or Reykjavik 
for which accurate long-term epidemiological data are available [5Hl2|. Remarkably, in some cases, more 
than four orders of magnitude of scaling have been found [7| . Other features of scale invariance have been 
reported for measles [7] and other infectious child diseases [13l[T4] , rabies and bovine tuberculosis [15] , or 
cholera [l6j . 

At a theoretical level, as pointed out by Rhodes and Anderson [9lfT0] , the lack of a characteristic scale 
in epidemic outbreaks is reminiscent of earthquakes and their associated (Guttenberg-Richter) power- 
law distribution. Actually, the presence of scale-invariance in measles, pertussis and others has been 
justified in |12| by exploiting the analogies between simple models for such epidemics and well-known 
(self-organized) earthquake models [T]-[3]. In what follows, we focus now on meningitis, for which a related 
self-organized mechanism has been recently proposed [T7H52]. 
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The bacteria responsible for meningitis, Neisseria meningiditis or meningococcus, is a human com- 
mensah it is typicahy harmless and it is present in up to one fourth of the human population |23j . 
Infection is transmitted through close contact with previously infected individuals. It is noteworthy that 
killing their hosts is a highly undesirable outcome for bacteria; therefore it makes sense that evolution 
selected for hardly harmful bacteria strains. Nevertheless, the meningococcus can accidentally mutate 
into a potentially dangerous strain, becoming highly damaging or even lethal for the host. This is an 
example of a more general type of "accidental pathogens" that innocuously cohabit with the host but 
that eventually -even if rarely- mutate causing symptomatic disease [23j . 

Aimed at modeling accidental pathogens and to shed some light on the reasons for the emergence of 
scale invariance in meningococcal epidemics, StoUenwerk, Jansen and coworkers proposed a simple and 
elegant mathematical model [T7H22j . The StoUenwerk- Jansen (SJ) model is a variant of the basic 
susceptible-infected-recovered-susceptible (SIRS) model. In the SIRS, individuals can be in any of the 
following states: "susceptible" (S), "infected" (I), or "recovered" (R) (also called "immune") [21H2S]- 
Perfectly mixed populations are usually considered, i.e. every individual is neighbor of any other (which 
is a mean-field assumption in the language of Statistical Physics or a "panmictic" one in the Ecology 
jargon). The additional key ingredient introduced in the SJ model is a second, potentially dangerous, 
strain of infected individuals labeled Y. Dangerous strains appear at a very small (mutation) rate at 
every contagion event. The dynamics of Y is almost identical to that of / except for the fact that at 
a certain rate e they can cause meningococcal disease of newly infected hosts and eventually kill them, 
y H> X (see below) 

By working out explicitly the analytical solution of this mean-field model as well as performing 
computer simulations, the authors above came to the counterintuitive conclusion that the smaller the 
value of e the larger the total amount of individuals killed on average in a given outbreak [T7l[T8]. This 
apparent paradox is easily resolved by realizing that the total number of individuals infected with Y grows 
upon decreasing e. Actually, it has been shown that, whilst for high pathogenicity e the distribution of 
the number of observed K-cases. s, is exponential, in the limit e — > it becomes power-law distributed, 
obeying 

Pis)^s~^-='^'^T{s/s,). (1) 

where is some scaling function and Sc a maximum characteristic scale controlled by 1/e [TTIITS]. 

Observe that the exponent r = 3/2 matches that of a critical branching process, for which the average 
population of infected individuals does not either grow or decay in time but stays constant |27| . Also, 
r = 3/2 coincides with the exponent for the distribution of first return times of an unbiased random- 
walk [28], which describes generically the scaling of avalanches in mean- field models. 

However, the value r = 3/2 predicted by the SJ model in its mean-field version, does not necessarily 
correspond to the best fit to empirical data [TICTTU] . To scrutinize the possible origin of such discrepancies, 
one would like to go beyond the mean-field/panmictic hypothesis by considering structured populations 
in which each individual has a finite local neighborhood. In what follows we shall study the SJ model 
on populations distributed on regular (Euclidean) lattices in dimensions d = 2 and d = 1. The study 
of more complex networks (as small- world networks or networks with communities), aimed at describing 
more realistically the net of social contacts, is left for a future work. 

As already pointed out [T71[T1], the SJ model can be straightforwardly made spatially explicit. Even 
if analytical or numerical calculations have not been performed, it was conjectured in |17H20| that in d- 
dimensions the SJ model should be in the directed percolation class, i.e. the broad and robust universality 
class characterizing phase transitions between an active and an absorbing state p9ti34j . In the present 
case, the absorbing state would be the F-free state obtained for sufficiently small F-infection rates, while 
the active or endemic phase would correspond to a non-vanishing density of Y's. In particular, if the 
critical behavior was indeed that of directed percolation, then r w 1.268 for two-dimensional populations 
(and T « 1.108 for epidemics propagating in one dimension) (see [35] and [36]). 

A priori, the prediction of directed percolation scaling is somehow suspicious if the model is indeed 
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self-organized; it has been shown that in general models of self-organized criticality, as sandpiles, rieepiles, 
earthquake models, etc. are not in the directed percolation class [57H41| . Furthermore, models of self- 
organized criticality lacking of any conservation law (as is the ease of the SJ model) have been shown not to 
be strictly critical, i.e. they are just approximately close to critical points [42] : instead the exact solution 
by StoUenwerk and Jansen proves that their model is exactly critical [ITlIT^. These considerations cast 
some doubts on the conjecture of the SJ model being a model of self-organized criticality |Tl[2| l37ti4Tj . 
Therefore, one is left with the following open questions: what is the SJ model true critical behavior?, 
what is the key ingredient why accidental pathogens - as described by the SJ model- originate scale 
invariant outbreaks?, does such an ingredient appear in other epidemics? 

The purpose of the present paper is to answer these questions by analyzing the SJ model in spatially 
extended systems. 



The spatially explicit SJ model is defined as follows. Each site of a d-dimensional lattice is in one of 
the following states S,I,Y,R, or X. The dynamics proceeds at any spatial location x and any time t 
according to the following one-site processes: 

• Spontaneous recovery of the benign strain: / — _R, at rate 7. 

• Spontaneous recovery of the dangerous strain: Y ^ R, at rate 7. 

• Loss of immunity: i? — ?► 5 at rate a. 

• Replacement or recovery of diseased: X S, at rate 
and two-site processes: 

• Infection with the benign strain: 5 + + /, at rate (3 — fi. 

• Infection with the dangerous strain: S + Y ^ Y + Y , at rate (3 — e. 

• Mutation: S + I ^ Y + I , at rate fi. 

• Disease: S + Y ^ X + Y, at rate e. 

Some comments on these reaction rules are in order. Accordingly to [51] it is assumed that being 
infected with one strain protects against co-infection with a second one. This assumption, which is 
supported by some empirical observations [21] . prevents transitions as I + Y ^ Y + Y from appearing. X 
(diseased/dead) individuals are immediately replaced by new susceptible ones, so that the total population 
size is kept fixed; i.e. 1^9 —> 00. A "back mutation" {S + Y ^ Y + I) reaction could be introduced, but 
for most purposes its rate is so small that it can be neglected. The dynamics is easily implemented in 
computer simulations by using the Gillespie's algorithm j43j or a variation of it in a rather standard way. 
In the well- mixed case writing the densities of the different species as S, /, Y, R, and X, one readily 
obtains the following mean-field or rate equations: 



Methods 
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aR + <j)X- PS{I + Y) 
iP - fi)SI - 7/ 
{13 - e)SY --/Y + fiSI 
"f{I + Y)~aR 
eSY - ipX 



(2) 
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satisiying the constraint S + I + Y + R + X^ 1. StoUenwerk and Jansen [TTIIT^ worked out the exact 
solution of this set of equations, concluding that it is critical in the limit e — > 0, and obtaining the explicit 
form of Eq.([T]). 

Instead, in spatially explicit models this mean-field approach breaks down: (i) the densities need to be 
replaced by spatio-temporal fields, as S{x,t), I{x,t), Y(x,t), etc; (ii) new (Laplacian) terms describing 
the map of nearest-neighbor contacts appear; (iii) fiuctuations become relevant and noise terms need 
to be added to account for demographic fluctuations. A full set of such stochastic Langevin equations 
can be derived from the microscopic dynamics by using standard procedures [28 [ I44 | |45 ] and numerically 
investigated [IB] (results not shown here). 

Results 

We now report on extensive numerical simulations for the SJ model. Wc have chosen the following 
parameter values: a = 0.3, (3 ~ 0.2 and 7 = 0.1, and have tried other different sets to confirm the 
robustness of the results. We consider the mutation rate fi to be sufficiently small, such that strains 
generated by consecutive mutations do not overlap; i.e. outbreaks finish well before a new mutation 
appears. For this reason we take as initial condition for any outbreak a state with a single mutant of the 
potentially dangerous strain Y and effectively fix = during outbreaks. 

Two dimensions 

We consider square lattices of linear size L — 32, 64, 128, 256, 512, and L = 1024, take a random initial 
condition with only S and / individuals, and run the dynamics with periodic boundary conditions, keeping 
/i = 0, until a steady state is reached. For instance, for L = 128, the steady state is characterized by 
< / >~ 0.38, < R >~ 0.12 and < 5* >~ 0.50 for the set of parameters above. 

Once the system sets into its steady state, we place a Y individual at the geometrical center of the 
lattice and study its spreading [3TH34] . To avoid finite size effects, spreading experiments are stopped 
once the K-strain touches the system boundary, and the described procedure is iterated. Depending 
on system size we ran up to 10^ independent realizations. As customarily done, we monitor: (i) the 
epidemics size distribution, analogous to Eq.([T]), for both Y and X; (ii) the average total number of Y 
as a function of time N{t); (iii) the surviving probability Psit) that the Y strain is still present in the 
system at time t, and (iv) the average square radius from the origin of y- infected individuals, R^{t). At 
criticality, these quantities are expected to scale algebraically as N{t) ~ t"^ , Ps{t) , and R^{t) ^ i^, 

while they should show exponential cut-offs in the sub-critical (or absorbing) phase. 

FiglT] shows the epidemic outbreak size distribution for different values of e and system size 128 * 128 
(size of y-infected outbreaks, sy, in the main plot, and size of X-infected outbreaks, Sx, in the inset). 
The probability distribution of avalanches sizes sx for X-infected sites is observed to inherit the statistics 
of y- infected ones. Both distributions can be well fitted by Eq.(IT]) where is a cut-off (exponential) 
function and Sc e"^/'^ determines the maximum size. The best fit gives r « 1.48(5) (fully compatible 
with 3/2 as shown in FiglT]) and 1/cr w 2 (not shown) both for P{sy) and P{sx)- From these plots we 
conclude that the system becomes critical in the limit of vanishing e, as occurs in mean-field. Observe 
that for small pathogeneicities, as e = 0.02, the system, even if sub-critical, exhibits scaling along more 
than three orders of magnitude. 

Both the mean and the variance of the above distributions are observed to diverge as power-laws in 
the limit e — ?• 0+. Actually, as shown in Figl2]the ratio of the variance ct^^ (resp. ag^) to the mean (sy) 
(resp. {sx}) diverges when e — > 0+ as (^syK^y) ~ <^sx/{^x) ~ e~^-^^^^ in the infinitely large system-size 
limit (otherwise, for finite L a size-induced cut-off appear as illustrated in Figj2]by comparing simulations 
for two different sizes). 
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Analogously, fixing e = wc measure the mean and variance as a function of system size; the ratio of 
the variance to the mean diverges very fast as L grows, I {sy) ~ (results not shown). Observe 

that, trivially, in this case P{sx) = S{sx), i-e. no death is produced. 

FiglSl shows results obtained for spreading quantities for various sizes and e = 0, i.e. at criticality. 
The best fits we obtain for the asymptotic behavior of these three magnitudes are 

N{t)^l, R'{t)^t (3) 

An effective power-law, with exponent slightly smaller than unity can be fit to our numerical results for 
Ps{t) (see FiglH]); however, such an effective value of the exponent can be seen to grow upon extending 
the maximum time, making the case for a logarithmic correction as described by Eq.([3]) and illustrated 
in FigUl Observe also that the reported values, 77 = 0, i5=l, z = l satisfy the hyper-scaling relation 
S + ri = dz/2 using d = 2 [3S]. 

Contrarily to a priori expectations and to a previous conjecture, these asymptotic laws have nothing 
to do with directed percolation values (which predicts pure power law behavior for the three spreading 
quantities [HHSl] with 77 « 0.229, S w 0.450, z w 1.132, r w 1.268; see [35] for a set of spreading and 
avalanche exponent numerical values in different universality classes). 

Instead, our results for spreading are in excellent agreement with the expectations for the two- 
dimensional voter-model [47H49] universality class (also called "compact directed percolation class" 
[311f33l [50]) [SUISSIET] (see below). Also, our results for the size distribution are in excellent agree- 
ment with the two-dimensional voter class expectations, r = 3/2 and 1/cr = 2. In particular, using these 
theoretical values, the variance to mean ratio should scale as L"* in good agreement with the numerical 
finding above. 

One dimension 

To further confirm the conclusion above, we have also performed studies of one-dimensional lattices, for 
which the expected behavior in the voter-model class is: 

7V(0^t°, Psit)-t-'/^ R\t)^t, (4) 

while the avalanche size exponent is t = 4/3 |35| . 

FiglS] shows the avalanche size distribution for various one dimensional sizes (512, 1024,2048,4096, 
and 8192) and e = 0. The measured r exponent, r = 1.33(3) is in excellent agreement with the voter- 
model class expectation. Fig|6] shows the result of spreading experiments (for L = 1024, 2048,4096, and 
8192). The best fits to the slopes, 77 « 0, 5 = 0.50(2), and z = 1.01(1) are in excellent agreement 
with the theoretical prediction (see dashed lines in Fig|6]). These results confirm that the SJ model is 
in the voter-class also in one spatial dimension (and exclude one-dimensional directed percolation values 
77 « 0.313, 5 « 0.159, z « 1.265, t « 1.108). 

Discussion 

Neutral theories date back to the sixties when Kimura introduced them in the context of population 
genetics [52j . Kimura assumed, as a null or neutral model, that each allele of a given gene (in haploid 
populations) is equally likely to enter the next generation, i.e. allele-type does not affect the prospects 
for survival or reproduction. Inspired in this, Hubbell proposed an analogous neutral theory of (forest) 
bio-diversity, in which the prospects of death and reproduction do not depend on the tree species |53j . 
Both theories lead to correct predictions and also to interpretation controversies (see [54] for a recent 
review) . 
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In its simplest spatially explicit version the neutral theory is as follows. Consider, for argument's sake, 
the neutral theory of bio-diversity with only two tree species. It can be formulated by considering the 
spatially explicit Moran process f55| or voter model [34,47 it9j . This is defined by two symmetrical species 
(up/down, right/left, black/white, etc.) fully occupying a d-dimensional lattice (or, more in general, an 
arbitrary network). At some rate a randomly chosen individual is removed and replaced by any of its 
nearest neighbors with homogeneous probability. This leads to a local tendency to create clusters of any 
of the two symmetrical species. Clusters of each of the two species will occupy different positions until 
eventually one of the two species will take over the whole (finite) space, leading to fiuctuation induced 
mono-dominance in any finite system-size. 

Such a type of coarsening process has been studied in depth; the asymptotic scaling of the spreading 
quantities, N{t), Ps{t), and R^{t) is analytically predicted to be given by Eq.Q. In particular, the 
logarithmic corrections for the surviving probability correspond to the fact d = 2 is the upper critical 
dimension of the voter universality class and are closely related to the marginality of the return time of 
two-dimensional random walkers |36[|47I - H^ . The key ingredients of this universality class turn out to be 
the symmetry between the two existing absorbing (mono-dominated) states [51 , 56, ST, and the absence 
of surface tension [51] . 

It is important to underline that, given the symmetry (neutrality) between the two species in the voter 
model, a given population of any of the species can either grow or decline with the same probability; i.e. 
there is no deterministic bias. Using the field theoretical jargon, the "mass" or "gap" term vanishes; all 
this is tantamount to the model being critical [55] . 

Indeed, for the voter model in mean-field, calling (j) the density in one of the two states (and 1 — (p 
the complementary density for the second species), then 

-0(1 - 0) + (1 - 0)0 = (5) 

where the first term represents the loss of an individual in the first state by contact with a neighboring 
in the second one, and the second represents the reverse process. Therefore, the average density of is 
constant all along the system evolution. Actually, in the voter-model there is no control parameter to be 
tuned; the model lies, by definition at a critical point: criticality is imposed by the neutral symmetry. 
Instead, introducing a bias or preference towards one of the species, 

4>^-{l- a)0(l - 0) + (1 + a){l - 0)0 = a{l - 0)0 (6) 

where the constant a > quantifies the degree of asymmetry, a non- vanishing linear or "mass" term is 
generated, bringing about an overall tendency for to grow or to decrease depending on the sign of a, 
i.e. deviating the system from criticality. 

Introducing spatial dependence and stochasticity, Eq.Q transforms into the Langevin equation for 
the voter class [51] 

0(x, t) - V20(x, t) + V<^(x, t){l - 0(x, t))7](x, t) (7) 

where 7y(x, t) is a Gaussian white noise (observe that the equation is symmetrical under the change 
<H- (1 — 0) and that there are two absorbing states at = and 0=1, respectively). 

Despite of the coincidence in the asymptotic scaling, let us underline that the SJ model is not a voter 
model. In the SJ model there are 5 different species and not just 2. In the case 6 = 0, however, / and Y 
are perfectly symmetrical; but as processes as F-t-I^-F-l-F or / do not exist, replacement 

of one species by the other occurs only if mediated by S particles. It is, therefore, only at sufficiently 
coarse grained scales that the dynamics behaves as the voter-model; microscopically the two dynamics 
differ significantly. 

Even if the SJ is not a voter model, the underlying reason for it to exhibit scale- invariance is that, in 
the limit e = and /i = the model is "neutral" (i.e. strains / and Y are perfectly symmetrical) and. 
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as a direct consequence, it is critical. More explicitly, calling Z the total density of infected sites which 
includes both / and Y then, at a mean field level, keeping fi ^ and e = 

S = aR- PSZ 
Z = PSZ - 7Z 

R = jZ-aR, (8) 

and the steady state is Zst = (a/3 — aj)/{af3 + f3j), Sst = 7//?, and Rst = Zst/y/a. Let us suppose, 
for argument sake, that Zst > 0; i.e. that there is a non-vanishing stationary density of infected sites 
(otherwise epidemics would just extinguish in finite time) and that at time i = is Y{t = 0) = Yq. The 
evolution of Y is controlled by 

Y = l3SstY --fY^O, (9) 

implying Y{t) = Yq, i.e. there is no bias for Y to grow or decay, or in other words, under these conditions 
the model is critical in what respects the Y variable. For example, in spreading experiments we fix 
Yq = 1/N (and Iq = Zgt — ^/N), implying that the total number of F-sites is 1 on average, in agreement 
with the numerical findings in Figl3] and FiglHl Switching on a non- vanishing e is equivalent in the 
voter model to introduce a bias, inducing deviations from criticality, i.e. exponential cut-offs for the size 
distribution and spreading quantities, as indeed observed in our numerical simulations for e ^ (see, for 
instance, Figl3]and Figl6]). 

It is noteworthy that in the perfectly symmetrical (neutral) case, e = 0, = 0, the model is somehow 
dull: two benign strains compete in a critical way, but there is no observable consequence of this for 
the population under study. The interesting behavior of the SJ model comes from slight deviations from 
criticality, this is, /i — > and e with e » fi; it is in this double limit that outbreaks leave behind an 
observable scale-invariant distribution of sick/dead individuals. 

In summary, the StoUenwerk-Janssen model provides us with the most parsimonious explanation 
for the appearance of scale-invariant epidemic outbreaks caused by accidental pathogens such as the 
meningococcus. Our main finding is that the system is critical in the limit of vanishing pathogenicity in 
all dimensions, and the reason for this is that the SJ is a neutral model, which turns out to be critical 
for the very same reasons as other neutral theories in Population Genetics and Ecology are critical. 
As a consequence of this, the critical behavior of the model is not described by directed percolation 
as previously conjectured, but instead by the voter-model universality class, representative of neutral 
theories. Understanding the theory of accidental pathogens as a neutral theory gives us new insight into 
the origin of generic scale invariance in epidemics in particular and in propagation phenomena in general. 

Similar ideas might apply to related problems as the spreading of computer viruses in the Internet or 
in the web of e-mail contacts: the largest overall damage is expected to occur for any type of spreading 
agents if their probability to cause damage is as small as possible. Being close to neutrality warrants 
success on the long term. 
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Figure 1. Main: Avalanche size distribution, P{s) for various values of e in a two dimensional lattice of size size 

128 * 128, and associated distribution for avalanche sizes of y-infected individuals, sy (main plot), and sx of X-infected 

ones (inset). 




Figure 2. Ratio of the variance (as) to the mean size (s) for Y (red) and for X (black), for and different values of and 
system sizes 128 * 128 (full symbols) and 256 * 256 (empty symbols). The ratio diverges in the double limit e — > 0, L — >■ oo. 
Inset; as the main plot, but in double logarithmic scale: as system-size increases the curves converge to straight lines, i.e. 
power-laws, but finite size effects are significant. 
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Figure 3. Log log plot of N{t), Ps{t), and R^{t), for spreading experiments in different two-dimensional systems (linear 
sizes L = 32, 64, 128, 256, 512, and 1024) as a function of time t. Dashed lines are a guide to the eye. 




Figure 4. Plot of tPs{t) as a function of log(t) for different system sizes (from L = 64 to L = 1024), illustrating the 
presence of linear logarithmic corrections for the surviving probability in the two-dimensional SJ model. 
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Figure 5. Avalanche size distribution, P(sy) for e = and various one-dimensional lattices, of size 

L = 512, 1024, 2048, 4096, and 8192, respectively. The straight line corresponds to the theoretical prediction for the 

onc-dimcnsional voter-model class. 




t 



Figure 6. Evolution of N{t), Ps{t), and R?{t) as a function of time t in log-log scale for spreading experiments 
performed on one-dimensional lattices of linear sizes 1024, 2048, 4096, and 8192. Parameter values are a = 0.3, 7 = 0. 
and /3 = 1.8. Results are averaged over 90000 independent realizations. 



